docs/notebooks/GEM-PRO - List of Gene IDs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GEM-PRO - List of Gene IDs\n",
    "\n",
    "This notebook gives an example of how to run the GEM-PRO pipeline with a **list of gene IDs**.\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "**Input:** List of gene IDs\n",
    "\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "**Output:** GEM-PRO model\n",
    "\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "import logging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import the GEM-PRO class\n",
    "from ssbio.pipeline.gempro import GEMPRO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Printing multiple outputs per cell\n",
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logging\n",
    "\n",
    "Set the logging level in `logger.setLevel(logging.<LEVEL_HERE>)` to specify how verbose you want the pipeline to be. Debug is most verbose.\n",
    "\n",
    "- `CRITICAL`\n",
    "     - Only really important messages shown\n",
    "- `ERROR`\n",
    "     - Major errors\n",
    "- `WARNING`\n",
    "     - Warnings that don't affect running of the pipeline\n",
    "- `INFO` (default)\n",
    "     - Info such as the number of structures mapped per gene\n",
    "- `DEBUG`\n",
    "     - Really detailed information that will print out a lot of stuff\n",
    "     \n",
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Warning:** \n",
    "`DEBUG` mode prints out a large amount of information, especially if you have a lot of genes. This may stall your notebook!\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create logger\n",
    "logger = logging.getLogger()\n",
    "logger.setLevel(logging.INFO)  # SET YOUR LOGGING LEVEL HERE #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Other logger stuff for Jupyter notebooks\n",
    "handler = logging.StreamHandler(sys.stderr)\n",
    "formatter = logging.Formatter('[%(asctime)s] [%(name)s] %(levelname)s: %(message)s', datefmt=\"%Y-%m-%d %H:%M\")\n",
    "handler.setFormatter(formatter)\n",
    "logger.handlers = [handler]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialization of the project\n",
    "\n",
    "Set these three things:\n",
    "\n",
    "- `ROOT_DIR`\n",
    "    - The directory where a folder named after your `PROJECT` will be created\n",
    "- `PROJECT`\n",
    "    - Your project name\n",
    "- `LIST_OF_GENES`\n",
    "    - Your list of gene IDs\n",
    "    \n",
    "A directory will be created in `ROOT_DIR` with your `PROJECT` name. The folders are organized like so:\n",
    "```\n",
    "    ROOT_DIR\n",
    "    └── PROJECT\n",
    "        ├── data  # General storage for pipeline outputs\n",
    "        ├── model  # SBML and GEM-PRO models are stored here\n",
    "        ├── genes  # Per gene information\n",
    "        │   ├── <gene_id1>  # Specific gene directory\n",
    "        │   │   └── protein\n",
    "        │   │       ├── sequences  # Protein sequence files, alignments, etc.\n",
    "        │   │       └── structures  # Protein structure files, calculations, etc.\n",
    "        │   └── <gene_id2>\n",
    "        │       └── protein\n",
    "        │           ├── sequences\n",
    "        │           └── structures\n",
    "        ├── reactions  # Per reaction information\n",
    "        │   └── <reaction_id1>  # Specific reaction directory\n",
    "        │       └── complex\n",
    "        │           └── structures  # Protein complex files\n",
    "        └── metabolites  # Per metabolite information\n",
    "            └── <metabolite_id1>  # Specific metabolite directory\n",
    "                └── chemical\n",
    "                    └── structures  # Metabolite 2D and 3D structure files\n",
    "                \n",
    "```\n",
    "\n",
    "<div class=\"alert alert-info\">**Note:** Methods for protein complexes and metabolites are still in development.</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# SET FOLDERS AND DATA HERE\n",
    "import tempfile\n",
    "ROOT_DIR = tempfile.gettempdir()\n",
    "\n",
    "PROJECT = 'genes_GP'\n",
    "LIST_OF_GENES = ['b0761', 'b0889', 'b0995', 'b1013', 'b1014', 'b1040', 'b1130', 'b1187', 'b1221', 'b1299']\n",
    "PDB_FILE_TYPE = 'mmtf'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: Creating GEM-PRO project directory in folder /tmp\n",
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: /tmp/genes_GP: GEM-PRO project location\n",
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: 10: number of genes\n"
     ]
    }
   ],
   "source": [
    "# Create the GEM-PRO project\n",
    "my_gempro = GEMPRO(gem_name=PROJECT, root_dir=ROOT_DIR, genes_list=LIST_OF_GENES, pdb_file_type=PDB_FILE_TYPE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mapping gene ID --> sequence\n",
    "\n",
    "First, we need to map these IDs to their protein sequences. There are 2 ID mapping services provided to do this - through **KEGG** or **UniProt**. The end goal is to map a UniProt ID to each ID, since there is a comprehensive mapping (and some useful APIs) between UniProt and the PDB.\n",
    "\n",
    "<p><div class=\"alert alert-info\">**Note:** You only need to map gene IDs using one service. However you can run both if some genes don't map in one service and do map in another!</div></p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Methods"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. automethod:: ssbio.pipeline.gempro.GEMPRO.kegg_mapping_and_metadata\n",
    "    :noindex:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7bcb759e902b4fa18c6bbe7d1e22db17",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: 10/10: number of genes mapped to KEGG\n",
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: Completed ID mapping --> KEGG. See the \"df_kegg_metadata\" attribute for a summary dataframe.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Missing KEGG mapping:  []\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>kegg</th>\n",
       "      <th>refseq</th>\n",
       "      <th>uniprot</th>\n",
       "      <th>pdbs</th>\n",
       "      <th>sequence_file</th>\n",
       "      <th>metadata_file</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>b0761</th>\n",
       "      <td>eco:b0761</td>\n",
       "      <td>NP_415282</td>\n",
       "      <td>P0A9G8</td>\n",
       "      <td>1B9M;1H9S;1B9N;1O7L;1H9R</td>\n",
       "      <td>eco-b0761.faa</td>\n",
       "      <td>eco-b0761.kegg</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b0889</th>\n",
       "      <td>eco:b0889</td>\n",
       "      <td>NP_415409</td>\n",
       "      <td>P0ACJ0</td>\n",
       "      <td>2GQQ;2L4A</td>\n",
       "      <td>eco-b0889.faa</td>\n",
       "      <td>eco-b0889.kegg</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b0995</th>\n",
       "      <td>eco:b0995</td>\n",
       "      <td>NP_415515</td>\n",
       "      <td>P38684</td>\n",
       "      <td>1ZGZ</td>\n",
       "      <td>eco-b0995.faa</td>\n",
       "      <td>eco-b0995.kegg</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b1013</th>\n",
       "      <td>eco:b1013</td>\n",
       "      <td>NP_415533</td>\n",
       "      <td>P0ACU2</td>\n",
       "      <td>4JYK;4XK4;4X1E;3LOC</td>\n",
       "      <td>eco-b1013.faa</td>\n",
       "      <td>eco-b1013.kegg</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b1014</th>\n",
       "      <td>eco:b1014</td>\n",
       "      <td>NP_415534</td>\n",
       "      <td>P09546</td>\n",
       "      <td>3E2Q;4JNZ;3E2R;4JNY;2GPE;4O8A;3E2S;2FZN;1TJ1;1...</td>\n",
       "      <td>eco-b1014.faa</td>\n",
       "      <td>eco-b1014.kegg</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            kegg     refseq uniprot  \\\n",
       "gene                                  \n",
       "b0761  eco:b0761  NP_415282  P0A9G8   \n",
       "b0889  eco:b0889  NP_415409  P0ACJ0   \n",
       "b0995  eco:b0995  NP_415515  P38684   \n",
       "b1013  eco:b1013  NP_415533  P0ACU2   \n",
       "b1014  eco:b1014  NP_415534  P09546   \n",
       "\n",
       "                                                    pdbs  sequence_file  \\\n",
       "gene                                                                      \n",
       "b0761                           1B9M;1H9S;1B9N;1O7L;1H9R  eco-b0761.faa   \n",
       "b0889                                          2GQQ;2L4A  eco-b0889.faa   \n",
       "b0995                                               1ZGZ  eco-b0995.faa   \n",
       "b1013                                4JYK;4XK4;4X1E;3LOC  eco-b1013.faa   \n",
       "b1014  3E2Q;4JNZ;3E2R;4JNY;2GPE;4O8A;3E2S;2FZN;1TJ1;1...  eco-b1014.faa   \n",
       "\n",
       "        metadata_file  \n",
       "gene                   \n",
       "b0761  eco-b0761.kegg  \n",
       "b0889  eco-b0889.kegg  \n",
       "b0995  eco-b0995.kegg  \n",
       "b1013  eco-b1013.kegg  \n",
       "b1014  eco-b1014.kegg  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# KEGG mapping of gene ids\n",
    "my_gempro.kegg_mapping_and_metadata(kegg_organism_code='eco')\n",
    "print('Missing KEGG mapping: ', my_gempro.missing_kegg_mapping)\n",
    "my_gempro.df_kegg_metadata.head()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. automethod:: ssbio.pipeline.gempro.GEMPRO.uniprot_mapping_and_metadata\n",
    "    :noindex:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[2018-02-05 18:12] [root] INFO: getUserAgent: Begin\n",
      "[2018-02-05 18:12] [root] INFO: getUserAgent: user_agent: EBI-Sample-Client/ (services.py; Python 3.6.3; Linux) Python-requests/2.18.4\n",
      "[2018-02-05 18:12] [root] INFO: getUserAgent: End\n",
      "[2018-02-05 18:12] [root] WARNING: status is not ok with Bad Request\n",
      "[2018-02-05 18:12] [root] WARNING: Results seems empty...returning empty dictionary.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b4dafc97239c48b085827c21552e0b80",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: 0/10: number of genes mapped to UniProt\n",
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: Completed ID mapping --> UniProt. See the \"df_uniprot_metadata\" attribute for a summary dataframe.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Missing UniProt mapping:  ['b1130', 'b0889', 'b1221', 'b0761', 'b0995', 'b1187', 'b1013', 'b1299', 'b1014', 'b1040']\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] WARNING: Empty dataframe\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>uniprot</th>\n",
       "      <th>reviewed</th>\n",
       "      <th>gene_name</th>\n",
       "      <th>kegg</th>\n",
       "      <th>refseq</th>\n",
       "      <th>num_pdbs</th>\n",
       "      <th>pdbs</th>\n",
       "      <th>ec_number</th>\n",
       "      <th>pfam</th>\n",
       "      <th>seq_len</th>\n",
       "      <th>description</th>\n",
       "      <th>entry_date</th>\n",
       "      <th>entry_version</th>\n",
       "      <th>seq_date</th>\n",
       "      <th>seq_version</th>\n",
       "      <th>sequence_file</th>\n",
       "      <th>metadata_file</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: [uniprot, reviewed, gene_name, kegg, refseq, num_pdbs, pdbs, ec_number, pfam, seq_len, description, entry_date, entry_version, seq_date, seq_version, sequence_file, metadata_file]\n",
       "Index: []"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# UniProt mapping\n",
    "my_gempro.uniprot_mapping_and_metadata(model_gene_source='ENSEMBLGENOME_ID')\n",
    "print('Missing UniProt mapping: ', my_gempro.missing_uniprot_mapping)\n",
    "my_gempro.df_uniprot_metadata.head()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. automethod:: ssbio.pipeline.gempro.GEMPRO.set_representative_sequence\n",
    "    :noindex:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f5babcb32874496590598b7e764208b8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: 10/10: number of genes with a representative sequence\n",
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: See the \"df_representative_sequences\" attribute for a summary dataframe.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Missing a representative sequence:  []\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>uniprot</th>\n",
       "      <th>kegg</th>\n",
       "      <th>pdbs</th>\n",
       "      <th>sequence_file</th>\n",
       "      <th>metadata_file</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>b0761</th>\n",
       "      <td>P0A9G8</td>\n",
       "      <td>eco:b0761</td>\n",
       "      <td>1B9M;1H9S;1B9N;1O7L;1H9R</td>\n",
       "      <td>eco-b0761.faa</td>\n",
       "      <td>eco-b0761.kegg</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b0889</th>\n",
       "      <td>P0ACJ0</td>\n",
       "      <td>eco:b0889</td>\n",
       "      <td>2GQQ;2L4A</td>\n",
       "      <td>eco-b0889.faa</td>\n",
       "      <td>eco-b0889.kegg</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b0995</th>\n",
       "      <td>P38684</td>\n",
       "      <td>eco:b0995</td>\n",
       "      <td>1ZGZ</td>\n",
       "      <td>eco-b0995.faa</td>\n",
       "      <td>eco-b0995.kegg</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b1013</th>\n",
       "      <td>P0ACU2</td>\n",
       "      <td>eco:b1013</td>\n",
       "      <td>4JYK;4XK4;4X1E;3LOC</td>\n",
       "      <td>eco-b1013.faa</td>\n",
       "      <td>eco-b1013.kegg</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b1014</th>\n",
       "      <td>P09546</td>\n",
       "      <td>eco:b1014</td>\n",
       "      <td>3E2Q;4JNZ;3E2R;4JNY;2GPE;4O8A;3E2S;2FZN;1TJ1;1...</td>\n",
       "      <td>eco-b1014.faa</td>\n",
       "      <td>eco-b1014.kegg</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      uniprot       kegg                                               pdbs  \\\n",
       "gene                                                                          \n",
       "b0761  P0A9G8  eco:b0761                           1B9M;1H9S;1B9N;1O7L;1H9R   \n",
       "b0889  P0ACJ0  eco:b0889                                          2GQQ;2L4A   \n",
       "b0995  P38684  eco:b0995                                               1ZGZ   \n",
       "b1013  P0ACU2  eco:b1013                                4JYK;4XK4;4X1E;3LOC   \n",
       "b1014  P09546  eco:b1014  3E2Q;4JNZ;3E2R;4JNY;2GPE;4O8A;3E2S;2FZN;1TJ1;1...   \n",
       "\n",
       "       sequence_file   metadata_file  \n",
       "gene                                  \n",
       "b0761  eco-b0761.faa  eco-b0761.kegg  \n",
       "b0889  eco-b0889.faa  eco-b0889.kegg  \n",
       "b0995  eco-b0995.faa  eco-b0995.kegg  \n",
       "b1013  eco-b1013.faa  eco-b1013.kegg  \n",
       "b1014  eco-b1014.faa  eco-b1014.kegg  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Set representative sequences\n",
    "my_gempro.set_representative_sequence()\n",
    "print('Missing a representative sequence: ', my_gempro.missing_representative_sequence)\n",
    "my_gempro.df_representative_sequences.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mapping representative sequence --> structure\n",
    "\n",
    "These are the ways to map sequence to structure:\n",
    "\n",
    "1. Use the UniProt ID and their automatic mappings to the PDB\n",
    "2. BLAST the sequence to the PDB\n",
    "3. Make homology models or \n",
    "4. Map to existing homology models\n",
    "\n",
    "You can only utilize option #1 to map to PDBs if there is a mapped UniProt ID set in the representative sequence. If not, you'll have to BLAST your sequence to the PDB or make a homology model. You can also run both for maximum coverage."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Methods"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. automethod:: ssbio.pipeline.gempro.GEMPRO.map_uniprot_to_pdb\n",
    "    :noindex:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: Mapping UniProt IDs --> PDB IDs...\n",
      "[2018-02-05 18:12] [root] INFO: getUserAgent: Begin\n",
      "[2018-02-05 18:12] [root] INFO: getUserAgent: user_agent: EBI-Sample-Client/ (services.py; Python 3.6.3; Linux) Python-requests/2.18.4\n",
      "[2018-02-05 18:12] [root] INFO: getUserAgent: End\n",
      "[2018-02-05 18:12] [root] WARNING: status is not ok with Bad Request\n",
      "[2018-02-05 18:12] [root] WARNING: Results seems empty...returning empty dictionary.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e7375542b2694874853b29384c4ecde5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: 0/10: number of genes with at least one experimental structure\n",
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: Completed UniProt --> best PDB mapping. See the \"df_pdb_ranking\" attribute for a summary dataframe.\n",
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] WARNING: Empty dataframe\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Empty DataFrame\n",
       "Columns: []\n",
       "Index: []"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Mapping using the PDBe best_structures service\n",
    "my_gempro.map_uniprot_to_pdb(seq_ident_cutoff=.3)\n",
    "my_gempro.df_pdb_ranking.head()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. automethod:: ssbio.pipeline.gempro.GEMPRO.blast_seqs_to_pdb\n",
    "    :noindex:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c7a4d6e929b74e029ac295f39cf39c31",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: Completed sequence --> PDB BLAST. See the \"df_pdb_blast\" attribute for a summary dataframe.\n",
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: 5: number of genes with additional structures added from BLAST\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pdb_id</th>\n",
       "      <th>pdb_chain_id</th>\n",
       "      <th>hit_score</th>\n",
       "      <th>hit_evalue</th>\n",
       "      <th>hit_percent_similar</th>\n",
       "      <th>hit_percent_ident</th>\n",
       "      <th>hit_num_ident</th>\n",
       "      <th>hit_num_similar</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>b0761</th>\n",
       "      <td>1b9n</td>\n",
       "      <td>B</td>\n",
       "      <td>1091.0</td>\n",
       "      <td>5.530720e-119</td>\n",
       "      <td>0.931298</td>\n",
       "      <td>0.931298</td>\n",
       "      <td>244</td>\n",
       "      <td>244</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b0761</th>\n",
       "      <td>1o7l</td>\n",
       "      <td>D</td>\n",
       "      <td>1089.0</td>\n",
       "      <td>1.096280e-118</td>\n",
       "      <td>0.931298</td>\n",
       "      <td>0.931298</td>\n",
       "      <td>244</td>\n",
       "      <td>244</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      pdb_id pdb_chain_id  hit_score     hit_evalue  hit_percent_similar  \\\n",
       "gene                                                                       \n",
       "b0761   1b9n            B     1091.0  5.530720e-119             0.931298   \n",
       "b0761   1o7l            D     1089.0  1.096280e-118             0.931298   \n",
       "\n",
       "       hit_percent_ident  hit_num_ident  hit_num_similar  \n",
       "gene                                                      \n",
       "b0761           0.931298            244              244  \n",
       "b0761           0.931298            244              244  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Mapping using BLAST\n",
    "my_gempro.blast_seqs_to_pdb(all_genes=True, seq_ident_cutoff=.9, evalue=0.00001)\n",
    "my_gempro.df_pdb_blast.head(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Downloading and ranking structures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Methods"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. automethod:: ssbio.pipeline.gempro.GEMPRO.pdb_downloader_and_metadata\n",
    "    :noindex:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Warning:** \n",
    "Downloading all PDBs takes a while, since they are also parsed for metadata. You can skip this step and just set representative structures below if you want to minimize the number of PDBs downloaded.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fd0319be0f394b809f7e68e6b1466bf4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: Updated PDB metadata dataframe. See the \"df_pdb_metadata\" attribute for a summary dataframe.\n",
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: Saved 15 structures total\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chemicals</th>\n",
       "      <th>description</th>\n",
       "      <th>experimental_method</th>\n",
       "      <th>mapped_chains</th>\n",
       "      <th>pdb_id</th>\n",
       "      <th>pdb_title</th>\n",
       "      <th>resolution</th>\n",
       "      <th>structure_file</th>\n",
       "      <th>taxonomy_name</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>b0761</th>\n",
       "      <td>NI</td>\n",
       "      <td>ModE (MOLYBDATE-DEPENDENT TRANSCRIPTIONAL REGU...</td>\n",
       "      <td>X-RAY DIFFRACTION</td>\n",
       "      <td>A;B</td>\n",
       "      <td>1b9m</td>\n",
       "      <td>REGULATOR FROM ESCHERICHIA COLI</td>\n",
       "      <td>1.75</td>\n",
       "      <td>1b9m.mmtf</td>\n",
       "      <td>Escherichia coli</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b0761</th>\n",
       "      <td>NI</td>\n",
       "      <td>MODE (MOLYBDATE DEPENDENT TRANSCRIPTIONAL REGU...</td>\n",
       "      <td>X-RAY DIFFRACTION</td>\n",
       "      <td>A;B</td>\n",
       "      <td>1b9n</td>\n",
       "      <td>REGULATOR FROM ESCHERICHIA COLI</td>\n",
       "      <td>2.09</td>\n",
       "      <td>1b9n.mmtf</td>\n",
       "      <td>Escherichia coli</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      chemicals                                        description  \\\n",
       "gene                                                                 \n",
       "b0761        NI  ModE (MOLYBDATE-DEPENDENT TRANSCRIPTIONAL REGU...   \n",
       "b0761        NI  MODE (MOLYBDATE DEPENDENT TRANSCRIPTIONAL REGU...   \n",
       "\n",
       "      experimental_method mapped_chains pdb_id  \\\n",
       "gene                                             \n",
       "b0761   X-RAY DIFFRACTION           A;B   1b9m   \n",
       "b0761   X-RAY DIFFRACTION           A;B   1b9n   \n",
       "\n",
       "                             pdb_title  resolution structure_file  \\\n",
       "gene                                                                \n",
       "b0761  REGULATOR FROM ESCHERICHIA COLI        1.75      1b9m.mmtf   \n",
       "b0761  REGULATOR FROM ESCHERICHIA COLI        2.09      1b9n.mmtf   \n",
       "\n",
       "          taxonomy_name  \n",
       "gene                     \n",
       "b0761  Escherichia coli  \n",
       "b0761  Escherichia coli  "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Download all mapped PDBs and gather the metadata\n",
    "my_gempro.pdb_downloader_and_metadata()\n",
    "my_gempro.df_pdb_metadata.head(2)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. automethod:: ssbio.pipeline.gempro.GEMPRO.set_representative_structure\n",
    "    :noindex:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4a7d5a92ffbb4f67843937ee1731b58f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: 5/10: number of genes with a representative structure\n",
      "[2018-02-05 18:12] [ssbio.pipeline.gempro] INFO: See the \"df_representative_structures\" attribute for a summary dataframe.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>id</th>\n",
       "      <th>is_experimental</th>\n",
       "      <th>file_type</th>\n",
       "      <th>structure_file</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>b0761</th>\n",
       "      <td>REP-1b9n</td>\n",
       "      <td>True</td>\n",
       "      <td>pdb</td>\n",
       "      <td>1b9n-A_clean.pdb</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b0889</th>\n",
       "      <td>REP-2gqq</td>\n",
       "      <td>True</td>\n",
       "      <td>pdb</td>\n",
       "      <td>2gqq-A_clean.pdb</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b1013</th>\n",
       "      <td>REP-4xk4</td>\n",
       "      <td>True</td>\n",
       "      <td>pdb</td>\n",
       "      <td>4xk4-A_clean.pdb</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b1187</th>\n",
       "      <td>REP-1h9t</td>\n",
       "      <td>True</td>\n",
       "      <td>pdb</td>\n",
       "      <td>1h9t-A_clean.pdb</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>b1221</th>\n",
       "      <td>REP-1rnl</td>\n",
       "      <td>True</td>\n",
       "      <td>pdb</td>\n",
       "      <td>1rnl-A_clean.pdb</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             id  is_experimental file_type    structure_file\n",
       "gene                                                        \n",
       "b0761  REP-1b9n             True       pdb  1b9n-A_clean.pdb\n",
       "b0889  REP-2gqq             True       pdb  2gqq-A_clean.pdb\n",
       "b1013  REP-4xk4             True       pdb  4xk4-A_clean.pdb\n",
       "b1187  REP-1h9t             True       pdb  1h9t-A_clean.pdb\n",
       "b1221  REP-1rnl             True       pdb  1rnl-A_clean.pdb"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Set representative structures\n",
    "my_gempro.set_representative_structure()\n",
    "my_gempro.df_representative_structures.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<StructProp REP-1h9t at 0x7f3880a275c0>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "{'_structure_dir': '/tmp/genes_GP/genes/b1187/b1187_protein/structures',\n",
       " 'chains': [<ChainProp A at 0x7f38830e1c18>],\n",
       " 'date': None,\n",
       " 'description': 'FATTY ACID METABOLISM REGULATOR PROTEIN',\n",
       " 'file_type': 'pdb',\n",
       " 'id': 'REP-1h9t',\n",
       " 'is_experimental': True,\n",
       " 'mapped_chains': ['A'],\n",
       " 'notes': {},\n",
       " 'original_structure_id': '1h9t',\n",
       " 'resolution': 3.25,\n",
       " 'structure_file': '1h9t-A_clean.pdb',\n",
       " 'taxonomy_name': 'ESCHERICHIA COLI'}"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Looking at the information saved within a gene\n",
    "my_gempro.genes.get_by_id('b1187').protein.representative_structure\n",
    "my_gempro.genes.get_by_id('b1187').protein.representative_structure.get_dict()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Saving your GEM-PRO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p><div class=\"alert alert-warning\">**Warning:** Saving is still experimental. For a full GEM-PRO with sequences & structures, depending on the number of genes, saving can take >5 minutes.</div></p>"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    ".. automethod:: ssbio.pipeline.gempro.GEMPRO.save_json\n",
    "    :noindex:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[2018-02-05 18:12] [root] WARNING: json-tricks: numpy scalar serialization is experimental and may work differently in future versions\n",
      "[2018-02-05 18:12] [ssbio.io] INFO: Saved <class 'ssbio.pipeline.gempro.GEMPRO'> (id: genes_GP) to /tmp/genes_GP/model/genes_GP.json\n"
     ]
    }
   ],
   "source": [
    "import os.path as op\n",
    "my_gempro.save_json(op.join(my_gempro.model_dir, '{}.json'.format(my_gempro.id)), compression=False)"
   ]
  }
 ],
 "metadata": {
  "_draft": {
   "nbviewer_url": "https://gist.github.com/342e44766d1dda637e975a4ef61dae19"
  },
  "anaconda-cloud": {},
  "celltoolbar": "Raw Cell Format",
  "gist": {
   "data": {
    "description": "GEM-PRO pipeline example for a list of gene IDs",
    "public": false
   },
   "id": "342e44766d1dda637e975a4ef61dae19"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  },
  "widgets": {
   "state": {
    "17ae537fadb545a796f7f00c4e01e5fb": {
     "views": [
      {
       "cell_index": 18
      }
     ]
    },
    "8388f9f0d2e9404fa72ac7fe92054d4a": {
     "views": [
      {
       "cell_index": 22
      }
     ]
    },
    "ddd32aebde664a5dbf0955fce77e77b9": {
     "views": [
      {
       "cell_index": 14
      }
     ]
    },
    "de8593f89bd14e57ae95734b4c36fc9a": {
     "views": [
      {
       "cell_index": 16
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}